home *** CD-ROM | disk | FTP | other *** search
/ Whiteline: Alpha / Whiteline Alpha.iso / progtool / modula2 / lpr / mathlib0.mod < prev    next >
Encoding:
Modula Implementation  |  1994-09-22  |  7.1 KB  |  329 lines
  1. IMPLEMENTATION MODULE MathLib0;
  2. FROM Terminal IMPORT WriteString,WriteLn;
  3.  
  4. (*      Methods and approximations are taken from 
  5.         HART et al, 
  6.         Computer Approximations, 
  7.         SIAM Series in Applied Mathematics,
  8.         John Wiley and Sons, 
  9.         New York, 1968
  10. *)
  11. TYPE    RealCard = RECORD
  12.                       CASE  : BOOLEAN OF
  13.                         FALSE: x : REAL |
  14.                         TRUE : y : CARDINAL
  15.                       END
  16.                     END;
  17.                     
  18. CONST   RedMax = 16;
  19.  
  20. VAR     x,y,z   : REAL;
  21.         NegSign : BOOLEAN;
  22.         shift   : RECORD 
  23.                     CASE : BOOLEAN OF 
  24.                       TRUE:  I : INTEGER | 
  25.                       FALSE: C : CARDINAL 
  26.                     END 
  27.                   END;
  28.         n       : CARDINAL;
  29.         Red     : ARRAY [1..RedMax] OF REAL;
  30.  
  31. PROCEDURE entier(x : REAL) : INTEGER;
  32.   BEGIN
  33.     NegSign:=x<0.0;
  34.     x:=ABS(x);
  35.     IF x>32767.0+FLOAT(ORD(NegSign)) THEN
  36.       WriteLn; 
  37.       WriteString('Error:  REAL too large in entier ');
  38.       WriteLn;
  39.       HALT;
  40.       RETURN 0
  41.     END;
  42.     IF NegSign THEN 
  43.       RETURN -(TRUNC(-x)+1) 
  44.     ELSE 
  45.       RETURN TRUNC(x) 
  46.     END;
  47.   END entier;
  48.  
  49. PROCEDURE real(x : INTEGER) : REAL;
  50.   BEGIN
  51.     IF x<0 THEN
  52.       RETURN -FLOAT(ABS(x))
  53.     ELSE
  54.       RETURN FLOAT(x)
  55.     END
  56.   END real;
  57.  
  58.  
  59. PROCEDURE ln(A : REAL) : REAL;
  60.   VAR b : RealCard;
  61. (* Hart function 2701 *)
  62.   BEGIN
  63.     IF A<=0.0 THEN
  64.       WriteLn;
  65.       WriteString('Error: ln of value<=0.0 ');
  66.       WriteLn;
  67.       HALT;
  68.       RETURN 0.0
  69.     END;
  70.     b.x:=A;
  71.     shift.I:=b.y DIV 80H;
  72.     DEC(shift.I,7EH);                 (*shift.I by power of 2*)
  73.     b.y:= 3F00H+(b.y MOD 80H);      (*b now in range 0.5..1.0*)
  74.     z:=b.x;
  75.     WHILE z<(1.0/1.4142135623) DO 
  76.       DEC(shift.I);
  77.       z:=2.0*z;
  78.     END;
  79.     z:=(z-1.0)/(z+1.0);
  80.     y:=z*z;
  81.     y:=z*((0.89554061525)*y-3.31355617479)/(y-1.65677797691);
  82.     y:= y+real(shift.I)*0.6931471805599453094172321;
  83.     RETURN y;
  84.   END ln;
  85.  
  86. PROCEDURE exp(A : REAL) : REAL;
  87. (* Hart function 1800*)
  88.   VAR b : RealCard;
  89.   BEGIN
  90.     IF A<0.0 THEN
  91.       b.x:=-A;
  92.       IF b.x>87.0 (*370.0*) THEN RETURN 0.0 END;
  93.       NegSign:=TRUE;
  94.     ELSE
  95.       b.x:=A;
  96.       IF A>87.0 (*370.0*) THEN
  97.         WriteLn;
  98.         WriteString('Error: value too large in exp');
  99.         WriteLn;
  100.         HALT;
  101.         RETURN 0.0
  102.       END;
  103.       NegSign:=FALSE;
  104.     END;
  105.     shift.I:=b.y DIV 80H;
  106.     DEC(shift.I,7BH);             (*shift.I by power of 2*)
  107.     IF shift.I>0 THEN
  108.       b.y:=3D80H+(b.y MOD 80H)  (*b now in range 0..1/32*)
  109.     ELSE
  110.       shift.I:=0
  111.     END;
  112.     y:=b.x*b.x;
  113.     z:=b.x*(0.08332586539095234*y+5.000089649149271528);
  114.     y:=1.0+(2.0*z)/(y+10.000179298301740601-z);
  115.     FOR n:=shift.I TO 1 BY -1 DO
  116.       y:=y*y;
  117.     END;    
  118.     IF NegSign THEN y:=1.0/y END;
  119.     RETURN y;
  120.   END exp;
  121.  
  122. PROCEDURE sqrt(A : REAL) : REAL;
  123.   VAR b   : RealCard;
  124.       exp : CARDINAL;
  125.   BEGIN
  126.     IF A=0.0 THEN RETURN A END;
  127.     IF A<0.0 THEN
  128.       WriteLn;
  129.       WriteString('Error: sqrt of value<0.0 ');
  130.       WriteLn;
  131.       HALT;
  132.       RETURN 0.0
  133.     END;
  134.     b.x:=A;
  135.     shift.I:=b.y DIV 80H;
  136.     DEC(shift.I,7EH);          (*shift.I by power of 2*)
  137.     IF ODD(shift.I) THEN 
  138.       exp:=3E80H;
  139.       INC(shift.I)
  140.     ELSE
  141.       exp:=3F00H
  142.     END;
  143.     b.y:=exp + b.y MOD 80H;  (*b.x now in [0..1/2]*)
  144.     x:=b.x;
  145.     y:=((((2.97530391*x)+2.02772463)*x+1.09542405E-1)*x+3.16235E-4)/
  146.       (((x+3.46399556)*x+6.41225367E-1)*x+9.408909E-3);
  147.         (*this is the first guess at a result*)
  148.     REPEAT
  149.       z:=0.5*(x/y-y);
  150.       y:=y+z;
  151.     UNTIL ABS(z)<1.0E-7 (*-15*);
  152.         (* now add size back in *)
  153.     IF shift.I<>0 THEN
  154.       b.x:=y;
  155.       INC(b.y,(shift.C DIV 2)*80H);
  156.       RETURN b.x
  157.     ELSE    
  158.       RETURN y;
  159.     END;
  160.   END sqrt;
  161.   
  162. PROCEDURE Reduce(VAR x :REAL);
  163.   BEGIN
  164.     IF x>Red[RedMax] THEN
  165.       WriteLn;
  166.       WriteString('Error: Too big argument in trig. function');
  167.       WriteLn;
  168.       HALT;
  169.       x:=0.0;
  170.       RETURN
  171.     END;
  172.     n:=1;
  173.     WHILE x>Red[n] DO INC(n) END;
  174.     IF n>1 THEN
  175.       FOR n:=n-1 TO 1 BY -1 DO
  176.         WHILE x>Red[n] DO x:=x-Red[n] END
  177.       END
  178.     END
  179.   END Reduce;
  180.  
  181. PROCEDURE sin(A : REAL) : REAL;
  182.   BEGIN
  183. (* Hart function no 3362 *)
  184.     IF A<0.0 THEN
  185.       NegSign:=TRUE;
  186.       x:=-A;
  187.     ELSE
  188.       NegSign:=FALSE;
  189.       x:=A;
  190.     END;
  191.     Reduce(x);
  192.     IF x>pi THEN
  193.        NegSign:=NOT NegSign;
  194.        x:=x-pi;
  195.     END;
  196.     IF x>(pi*0.5) THEN
  197.       x:=pi-x;
  198.     END;
  199.     x:=(2.0/pi)*x;
  200.     y:=x*x;
  201.     y:= x*((10.42302058487*y-173.89497132272)*y+529.30152776255)/
  202.         ((y+27.86575519992)*y+336.96381989527);
  203.     IF NegSign THEN 
  204.       RETURN -y
  205.     ELSE 
  206.       RETURN y 
  207.     END;
  208.   END sin;
  209.  
  210. PROCEDURE cos(A : REAL) : REAL;
  211.   BEGIN
  212.     RETURN sin(A+pi/2.0);   (*can do better - eg Hart 3502*)
  213.   END cos;
  214.  
  215. PROCEDURE tan(A : REAL) : REAL;
  216.   PROCEDURE Hart4282(A : REAL) : REAL;
  217.     BEGIN
  218.       x:=4.0/pi*A;
  219.       y:=x*x;
  220.       RETURN x*(-12.55329742424*y+212.42445758263)/
  221.              ((y-71.59606050466)*y+270.46722349399);
  222.     END Hart4282;
  223.   BEGIN
  224.     IF A<0.0 THEN
  225.       x:=-A;
  226.       NegSign:=TRUE
  227.     ELSE
  228.       x:=A;
  229.       NegSign:=FALSE
  230.     END;
  231.     Reduce(x);
  232.     IF x>pi THEN x:=x-pi END;
  233.     IF x>(pi/2.0) THEN
  234.       x:=pi-x;
  235.       NegSign:=NOT NegSign;
  236.     END;
  237.     IF x>(pi/4.0) THEN 
  238.       z:=x;
  239.       x:=(pi/2.0-x);
  240.       IF x<1.0E-19 THEN 
  241.         y:=1.0E19 
  242.       ELSE    
  243.         z:=Hart4282(x);
  244.         y:=1.0/z;
  245.       END;
  246.     ELSE 
  247.       y:=Hart4282(x) 
  248.     END;
  249.     IF NegSign THEN 
  250.       RETURN -y
  251.     ELSE 
  252.       RETURN y 
  253.     END;
  254.   END tan;
  255.  
  256. PROCEDURE arctan(A : REAL) : REAL;
  257.   VAR NegSign:BOOLEAN;
  258.   PROCEDURE Hart5071(x : REAL) : REAL;
  259.     BEGIN
  260.       y:=x*x;
  261.       RETURN x*(6.3691887127*y+12.65998609915)/
  262.            ((y+10.5891113168)*y+12.65998646243);
  263.     END Hart5071;
  264.   BEGIN
  265.     IF A<0.0 THEN
  266.       NegSign:=TRUE;
  267.       A:=-A;
  268.     ELSE
  269.       NegSign:=FALSE;
  270.     END;
  271.     IF A>1.0E-4 THEN 
  272.        (* Reduce range to be in [0..tan(pi/8)]*)
  273.       IF A=1.0 THEN 
  274.         A:=pi/4.0
  275.       ELSIF A>1.0 THEN
  276.         A:=pi/2.0-arctan(1.0/A)
  277.       ELSIF A>=0.4142135623730950488 (*tan(pi/8)*) THEN
  278.         A:=pi/4.0-Hart5071(2.0/(1.0+A)-1.0)
  279.       ELSE
  280.         A:=Hart5071(A);
  281.       END;
  282.     END;
  283.     IF NegSign THEN 
  284.       RETURN -A
  285.     ELSE 
  286.       RETURN A
  287.     END;
  288.   END arctan;
  289.  
  290. PROCEDURE arctan2(Y, X : REAL) : REAL;
  291.   VAR Quadrant  :CARDINAL; 
  292.   BEGIN
  293.     IF (X>=0.0) THEN
  294.       IF Y>=0.0 THEN 
  295.         Quadrant:=1
  296.       ELSE
  297.         Y:=-Y;
  298.         Quadrant:=4;
  299.       END
  300.     ELSE
  301.       X:=-X;
  302.       IF Y>=0.0 THEN 
  303.         Quadrant:=2
  304.       ELSE
  305.         Y:=-Y;
  306.         Quadrant:=3;
  307.       END;
  308.     END;
  309.     IF X<=1.0E-19*Y THEN 
  310.       x:=pi/2.0
  311.     ELSE x:=arctan(Y/X) END;
  312.     CASE Quadrant OF
  313.       1: RETURN x    |
  314.       2: RETURN pi-x |
  315.       3: RETURN pi+x |
  316.       4: RETURN 2.0*pi-x
  317.     END;
  318.   END arctan2;
  319.  
  320. BEGIN  (* initialisation *)
  321.  
  322.   Red[1]:=2.0*pi;
  323.   FOR n:=2 TO RedMax DO
  324.     Red[n]:=3.0*Red[n-1]
  325.   END 
  326.  
  327. END MathLib0.
  328.  
  329.